Ειδικά Μαθήματα Βοτανικής

1 Βασικές προϋποθέσεις

Εννοείται ότι πριν ξεκινήσουμε να κάνουμε το οτιδήποτε, έχουμε δημιουργήσει ένα νέο project στο R-Studio, το όνομα του οποίου - για λόγους τους οποίους θα καταλάβετε αργότερα - ΔΕΝ πρέπει να περιέχει:
1. ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)

Το ίδιο ισχύει ΚΑΙ για το file path του εν λόγω φακέλου (π.χ., όχι ‘E:/Ειδικά Μαθήματα Βοτανικής/SDM’, αλλά ’E:/Eidika_Mathimata_Botanikis/SDM).

Ένα άλλο κρίσιμο σημείο είναι το εξής: θα χρειαστεί να εγκαταστήσετε τα πακέτα που πρόκειται να χρησιμοποιήσετε σε αυτό το tutorial1.

2 Εγκατάσταση και φόρτωση βιβλιοθηκών

Όπως κάθε φορά, θα χρειαστεί να εγκαταστήσουμε και να φορτώσουμε τις απαραίτητες βιβλιοθήκες.

## ===========================================================================##
## Install the main packages
## ===========================================================================##
install.packages(c("tidyverse", "rgbif", "raster", "data.table", "rgeos", "rgdal", 
    "gridExtra", "dismo", "biomod2", "sf", "usdm", "speciesgeocodeR", "pacman", 
    "spThin", "CoordinateCleaner", "ggspatial", "scales", "scrubr", "mapr", 
    "rasterVis", "ade4", "viridis", "biogeo"), dependencies = T)
## ===========================================================================##

Ας τις φορτώσουμε:

## ===========================================================================##
## Load them
## ===========================================================================##
pacman::p_load(tidyverse, rgbif, raster, usdm, rgeos, rgdal, usdm, biomod2, 
    gridExtra, countrycode, devtools, sf, spThin, biogeo, speciesgeocodeR, pacman, 
    dismo, CoordinateCleaner, ggspatial, scales, scrubr, mapr, rasterVis, data.table, 
    ade4, viridis)
## ===========================================================================##

3 Διανυσματικά δεδομένα

3.1 Σύνορα χωρών

Πρώτα θα κατεβάσουμε τις γεωγραφικές μεταβλητές για την χώρα που μας ενδιαφέρει (η Ελλάδα στην προκειμένη περίπτωση). Θα πρέπει να εισάγουμε (στην εντολή country στον κώδικα που παρατίθεται πιο κάτω) τον 3-ψήφιο ISO κωδικό για την χώρα μας που τον βρίσκουμε από αυτόν τον ιστότοπο ή με την εντολή ISO_countries <- getData("ISO3") %>% as.data.frame2. Εάν θέλουμε να περιορίσουμε την περιοχή μελέτης σε μια μικρότερη περιοχή, τότε θα πρέπει να χρησιμοποιήσουμε ένα άλλο λογισμικό, το Qgis, ώστε να “κόψουμε” την περιοχή που μας ενδιαφέρει. Προς το παρόν, δεν θα χρειαστεί να το κάνουμε αυτό, αλλά θα περιοριστούμε να κατεβάσουμε τα πλήρη δεδομένα για την Ελλάδα: .

## ===========================================================================##
## Download the data for Greece
## ===========================================================================##
## There are three levels available, you are free to explore them
Greece <- getData("GADM", country = "GRC", level = 3)
## ===========================================================================##

Μπορούμε να αναπαραστήσουμε γραφικά τα δεδομένα αυτά.

## ===========================================================================##
## Plot Greece
## ===========================================================================##
plot(Greece)

Είναι αρκετά εύκολο στην προκειμένη περίπτωση να περιορίσουμε την έκταση του πολυγώνου στην Κρήτη, με τις ακόλουθες εντολές:

## ===========================================================================##
## First we need to see what is inside the object we created
## ===========================================================================##
Greece
## class       : SpatialPolygonsDataFrame 
## features    : 326 
## extent      : 19.37236, 29.6457, 34.80069, 41.74801  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 16
## names       : GID_0, NAME_0,   GID_1,                      NAME_1,                                                                                        NL_NAME_1,     GID_2,         NAME_2,                        NL_NAME_2,       GID_3,   NAME_3,          VARNAME_3,                                                     NL_NAME_3, TYPE_3,    ENGTYPE_3, CC_3, ... 
## min values  :   GRC, Greece, GRC.1_1,                      Aegean,                                                                 <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1_1,          Athos, <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1.1_1,   Abdera, Abelokipi-Menemeni,                              <U+0386><U+03B8><U+03C9><U+03C2>,  Dímos, Municipality,   NA, ... 
## max values  :   GRC, Greece, GRC.8_1, Thessaly and Central Greece, Tessa<U+03BB><U+03AF>a <U+03BA>a<U+03B9> Ste<U+03C1>e<U+03AC> <U+0395><U+03BB><U+03BB><U+03AC>da, GRC.8.2_1, West Macedonia,           Tessa<U+03BB><U+03AF>a, GRC.8.2.9_1, Zografou,           Zografos, Tessa<U+03BB><U+03BF><U+03BD><U+03AF><U+03BA><U+03B7><U+03C2>,  Dímos, Municipality,   NA, ...
names(Greece)
##  [1] "GID_0"     "NAME_0"    "GID_1"     "NAME_1"    "NL_NAME_1"
##  [6] "GID_2"     "NAME_2"    "NL_NAME_2" "GID_3"     "NAME_3"   
## [11] "VARNAME_3" "NL_NAME_3" "TYPE_3"    "ENGTYPE_3" "CC_3"     
## [16] "HASC_3"
Crete <- subset(Greece, NAME_1 == "Crete")
## ===========================================================================##
## Plot Crete
## ===========================================================================##
plot(Crete)

3.2 GBIF

Στο σημείο αυτό, είμαστε πλέον στην ευχάριστη θέση να μπορούμε να χρησιμοποιήσουμε τα δεδομένα θέσης για τα είδη τα οποία μας ενδιαφέρουν. Υπάρχουν δύο περιπτώσεις, όσον αφορά την απόκτηση τέτοιων δεδομένων:
- να τα έχουμε αποκτήσει πρωτογενώς (δηλαδή να έχουμε εργαστεί στο πεδίο και να έχουμε εντοπίσεις τις θέσεις εμφάνισης του εκάστοτε είδους και των υποπληθυσμών του)
- να τα έχουμε αποκτήσει δευτερογενώς (δηλαδή είτε να έχουμε χρησιμοποιήσει τις διαδικτυακές βάσεις ψηφιοποιημένων - και ενίοτε γεωαναφερμένων - δεδομένων GBIF και BIEN, είτε να έχουμε ψηφιοποιήσει χάρτες εμφάνισης θέσης3)

Το ιδανικό βεβαίως είναι να έχουμε αποκτήσει πρωτογενώς τα στοιχεία αυτά, καθότι στην περίπτωση αυτή, μπορούμε να είμαστε σίγουροι για τις γεωγραφικές συντεταγμένες των υποπληθυσμών του προς μελέτη είδους. Στην πλειονότητα των περιπτώσεων όμως, δεν έχουμε πρόσβαση σε τέτοια δεδομένα, οπότε αναγκαστικά καταφεύγουμε στην δεύτερη λύση και δη στις διαδικτυακές βάσεις δεδομένων. Ευτυχώς για εμάς, υπάρχουν κάποια πακέτα4 στην R, όπως το rgbif και το rbien, τα οποία αυτοματοποιούν την διαδικασία μεταφόρτωσης δεδομένων από τις διαδικτυακές αυτές βάσεις5.

Εμείς θα εργαστούμε με το πακέτο rgbif προς το παρόν.

Σήμερα – όπως και στην προηγούμενη εργαστηριακή άσκηση – θα ασχοληθούμε ένα εκ των έξι μονοτυπικών γενών της ελληνικής χλωρίδας, την Petromarula pinnata.

Fig 1. The endemic species Petromarula pinnata from Crete.

Fig 1. The endemic species Petromarula pinnata from Crete.

Πρώτα χρειάζεται να βεβαιωθούμε ότι το προς μελέτη είδος, υπάρχει όντως εντός της διαδικτυακής βάσης την οποία επιθυμούμε να χρησιμοποιήσουμε, προκειμένου να εξάγουμε δεδομένα. Η εντολή name_suggest() ερευνά όλα τα είδη (taxa καλύτερα) τα οποία μοιάζουν με το είδος το οποίο ψάχνουμε και κρατά μόνο το όνομα των taxa το οποίο ταιριάζει ακριβώς με αυτό το οποίο μας ενδιαφέρει6.

spp_petromarula <- name_suggest(q = "Petromarula pinnata", rank = "species", 
    limit = 10000)
(spp_petromarula <- spp_petromarula[grepl("^Petromarula pinnata$", spp_petromarula$canonicalName), 
    ])
key canonicalName rank
3167989 Petromarula pinnata SPECIES

Στο βήμα αυτό αναγνωρίσαμε τον μοναδικό αριθμό αναφοράς (key) του προς μελέτη είδους και μόνο χρησιμοποιώντας τον αριθμό αυτό μπορούμε να είμαστε βέβαιοι ότι τα δεδομένα θέσης τα οποία θα λάβουμε από την διαδικτυακή αυτή βάση δεδομένων θα αναφέρονται στο είδος το οποίο μας ενδιαφέρει.

Στη συνέχεια, θα δούμε πόσα γεωαναφερμένα δεδομένα θέσης έχουμε στην διάθεση μας.

3.3 Μεταφόρτωση των γεωαναφερμένων δεδομένων θέσης

## ===========================================================================##
## Download data
## ===========================================================================##
sp_occ <- occ_search(taxonKey = spp_petromarula$key, country = "GR", hasCoordinate = T, 
    limit = 1000, return = "data")
## ===========================================================================##


## ===========================================================================##
## To prevent any problem with the pathway it is a good practice to remove
## blank spaces from species names
## ===========================================================================##
sp_occ$name <- sub(" ", "_", sp_occ$name)
## ===========================================================================##


## ===========================================================================##
## Total number of occurrences
## ===========================================================================##
sort(table(sp_occ$name), decreasing = T)
## Petromarula_pinnata A.DC. 
##                        77

Θα πρέπει να τονίσουμε το εξής: Εφ’όσον δεν έχουμε συλλέξει πρωτογενώς τα δεδομένα μας, θα χρειαστεί να ελέγξουμε εάν οι θέσεις εμφάνισης που μεταφορτώσαμε από τις διαδικτυακές βάσεις δεδομένων, συμφωνούν με τα γνωστά όρια εξάπλωσης του είδους που μας ενδιαφέρει. Εάν όντως υπάρχουν κάποια σημεία τα οποία βρίσκονται εκτός της περιοχής εξάπλωσης, θα χρειαστεί να τα αφαιρέσουμε από την ανάλυση.

## ===========================================================================##
## Clean the coordinates
## ===========================================================================##
flags_gbif <- clean_coordinates(x = sp_occ, lon = "decimalLongitude", lat = "decimalLatitude", 
    countries = "countryCode", species = "species", tests = c("capitals", "centroids", 
        "equal", "gbif", "zeros", "seas"), seas_ref = Crete)
## ===========================================================================##


## ===========================================================================##
## Keep only those that are OK
## ===========================================================================##
flags_gbif <- flags_gbif %>% mutate(Data = ifelse(.summary == "TRUE", "Clean", 
    "Flagged"), NAME_1 = "Crete", Dataset = "GBIF", longitude = decimalLongitude, 
    latitude = decimalLatitude, scientific_name = "Petromarula pinnata") %>% 
    dplyr::filter(str_detect(Data, "Cle")) %>% dplyr::select(scientific_name, 
    longitude, latitude, Dataset)
## ===========================================================================##


## ===========================================================================##
## Create a duplicate for safekeeping
## ===========================================================================##
flags_gbif_dpl <- flags_gbif
## ===========================================================================##

Ας φτιάξουμε έναν φάκελο ώστε να αποθηκεύσουμε τα δεδομένα μας και ας τα σώσουμε, προκειμένου να μπορέσουμε να τα φορτώσουμε κάποια άλλη στιγμή:

dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Petromarula pinnata", recursive = TRUE, showWarnings = FALSE)

saveRDS(flags_gbif, "./Species/ Petromarula pinnata.rds")

4 Μεταφόρτωση περιβαλλοντικών μεταβλητών

Τα κλιματικά δεδομένα μας μπορούμε να τα κατεβάσουμε από το WorldClim, ενώ τα υψομετρικά δεδομένα από εδώ και τα δεδομένα σχετικά με την ξηρασία και την εξατμισιοδιαπνοή από εδώ. - Εάν φυσικά δουλεύουμε για μια συγκεκριμένη εργασία (πτυχιακή, μεταπτυχιακή, κτλ) και όχι στα πλαίσια ενός προπτυχιακού μαθήματος επιλογής. Προς το παρόν, θα χρησιμοποιήσουμε κάποιες εντολές για να κατεβάσουμε τα δεδομένα που θέλουμε μέσω της R.

Όπως μπορείτε να διαπιστώσετε, υπάρχουν πάρα πολλά κλιματικά μοντέλα (επί παραδείγματι, CCSM4, HadGEM-2, κτλ) και τέσσερα (4) διαφορετικά κλιματικά σενάρια που έχουν να κάνουν με τη συγκέντρωση των αερίων του θερμοκηπίου στην ατμόσφαιρα εντός του 21ου αιώνα { ενδεικτική βιβλιογραφία ή άλλη ενδεικτική βιβλιογραφία ή ακόμα μια ενδεικτική βιβλιογραφία ή τελευταία ενδεικτική βιβλιογραφία}7.

Συνεπώς, δεν μπορούμε να χρησιμοποιήσουμε μόνο ένα κλιματικό μοντέλο, πόσω μάλλον ένα κλιματικό σενάριο. Στο σημείο αυτό, είναι εύλογο να αρχίσετε να πανικοβάλεστε8. Ευτυχώς όμως, δεν χρειάζεται και ούτε είναι επιστημονικά σωστό, να χρησιμοποιήσουμε όλα τα κλιματικά μοντέλα. Αναλόγως της περιοχής στην οποία εργαζόμαστε, κάποια κλιματικά μοντέλα αναπαριστούν καλύτερα τις περιβαλλοντικές συνθήκες σε σχέση με κάποια άλλα (McSweeney et al., 2015). Καθώς θα διαπιστώσετε (όταν με το καλό διαβάσετε την προαναφερθείσα εργασία), ακόμα και έτσι, είναι πολλά τα κλιματικά μοντέλα για την περιοχή μας. Όμως, μόνο 3 από αυτά έχουν δεδομένα και για τα 4 κλιματικά σενάρια (BCC, HadGEM, CCSM4). Εμείς χάριν συντομίας, θα εργαστούμε όμως μόνο με ένα (1) κλιματικό μοντέλο και ενα (1) κλιματικό σενάριο.

## ============================================================================##
## First, let's download current climate data for the whole planet
## ============================================================================##
current_climate <- getData("worldclim", var = "bio", res = 2.5)
## ============================================================================##


## ============================================================================##
## Then, let's download future climate data
## ============================================================================##
future_climate <- getData("CMIP5", var = "bio", res = 2.5, rcp = 85, model = "CC", 
    year = 70)

## ============================================================================##

Ας αλλάξουμε τα ονόματα των περιβαλλοντικών μεταβλητών.

names(current_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

names(future_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

Καθώς τα κλιματικά μας δεδομένα όπως διαπιστώσατε έχουν εξαιρετικά μεγάλο εύρος τιμών, θα χρειαστεί να τα μετατρέψουμε σε oC όσον αφορά την θερμοκρασία (bio_1). Εδώ μπορείτε να διαπιστώσετε το γιατί9. Αυτό γίνεται με την ακόλουθη εντολή:

gain(current_climate$bio_1) = 0.1

gain(future_climate$bio_1) = 0.1

Στη συνέχεια, θα κατεβάσουμε τα υψομετρικά δεδομένα για την χώρα που μας ενδιαφέρει:

altitude <- getData("alt", country = "GRC", mask = TRUE)

4.1 Κοπή των περιβαλλοντικών δεδομένων στην έκταση της Ελλάδας

Στη συνέχεια, θα ‘κόψουμε’ τις περιβαλλοντικές μας μεταβλητές στο σχήμα της Ελλάδας και ακολούθως θα δημιουργήσουμε ένα νέο αρχείο που θα περιέχει τις μεταβλητές αυτές και την περιοχή μελέτης μας:

## ============================================================================##
## First for Greece
## ============================================================================##
current_climate_crop <- crop(current_climate, Greece, snap = "in") %>% mask(., 
    Greece)

altitude <- crop(altitude, Greece, snap = "in") %>% mask(., Greece) %>% 
## We are forcing the extent of altitude to be exactly the same as that of
## the climate data
resample(., current_climate_crop, "bilinear")

future_climate <- crop(future_climate, Greece, snap = "in") %>% mask(., Greece)

## ============================================================================##


## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude)
future_climate_crop <- stack(future_climate, altitude)
## ============================================================================##


## ============================================================================##
## Then for Crete
## ============================================================================##
current_climate_Crete <- crop(current_climate_crop, Crete, snap = "in") %>% 
    mask(., Crete)

future_climate_Crete <- crop(future_climate_crop, Crete, snap = "in") %>% mask(., 
    Crete)

## ============================================================================##

4.2 Έλεγχος συγγραμικότητας και Variance Inflation Factor (VIF)

Όπως θα γνωρίζετε, ΔΕΝ μπορούμε να χρησιμοποιήσουμε σε μια ανάλυση δύο μεταβλητές, οι οποίες εμφανίζουν συντελεστή συσχέτισης \(> 0.7\)10, καθότι στην προκειμένη περίπτωση οι δύο αυτές μεταβλητές είναι συγγραμικές και ως εκ τούτου εάν τις συμπεριλάβουμε και τις δύο στην ανάλυση μας, τότε δεν θα γνωρίζουμε σε ποιά εκ των δύο θα οφείλεται το αποτέλεσμα της ανάλυσης μας.

Προκειμένου να ξεπεράσουμε αυτόν τον σκόπελο, θα τρέξουμε την ακόλουθη εντολή11:
.

uncorrelated_vars <- vifcor(current_climate_Crete %>% as.data.frame(), th = 0.7)
uncorrelated_vars <- uncorrelated_vars@results$Variables[uncorrelated_vars@results$VIF <= 
    5]
uncorrelated_vars <- droplevels(uncorrelated_vars)

current_climate_Crete <- stack(subset(current_climate_Crete, levels(uncorrelated_vars)))

future_climate_Crete <- stack(subset(future_climate_Crete, levels(uncorrelated_vars)))

5 Δεδομένα θέσης και περιβαλλοντικές μεταβλητές

Πλέον μπορούμε να συνδέσουμε τις γεωαναφερμένες θέσεις εμφάνισης του είδους που μας ενδιαφέρει με δεδομένα ψηφιδοπλέγματος. Στο σημείο αυτό, μπορούμε να εξάγουμε τα περιβαλλοντικά δεδομένα για κάθε μια εκ των θέσεων εμφάνισης του είδους που μας ενδιαφέρει με μια μόνο εντολή, αναδεικνύοντας κατ’αυτόν τον τρόπο το πλεονέκτημα του να έχουμε ενώσει προηγουμένως τα raster αρχεία μας σε ένα RasterStack.

coordinates(flags_gbif) <- c("longitude", "latitude")

pts.clim <- raster::extract(current_climate_Crete, flags_gbif, method = "bilinear")

sp_occ_clim <- data.frame(cbind(coordinates(flags_gbif), pts.clim, flags_gbif@data))

Και τώρα, ήρθε η στιγμή να οπτικοποιήσουμε τα αποτελέσματα μας και δη, την παρουσία του είδους μας στην περιοχή μελέτης:

map.ext <- extent(Crete)

plot(current_climate_crop[[20]], col = terrain.colors(100), ext = map.ext)

plot(flags_gbif, pch = 16, cex = 0.5, add = T)

plot(Crete, add = T)

Στο σημείο αυτό, θα ήθελα να σας τονίσω ότι ορισμένες φορές, πέρα από τους ελέγχους που διενεργήσαμε έως τώρα, ίσως χρειαστεί να προχωρήσουμε σε μια αραίωση (thinning) των θέσεων εμφάνισης. Αυτό μπορεί να είναι απαραίτητο όταν έχουμε πάρα πολλά σημεία τα οποία είναι αρκετά κοντά το ένα στο άλλο12. Σε κάθε περίπτωση, το επιστημονικά ορθό είναι να υπάρχει μια (1) θέση εμφάνισης ανά km2. Ένα πακέτο που μπορεί να μας βοηθήσει είναι το spThin13.

## Min. distance: 1000 m Number of reps: 1000
spoccs <- flags_gbif_dpl %>% as.data.frame() %>% dplyr::rename(x = longitude, 
    y = latitude, Species = scientific_name) %>% addmainfields(species = "Species")

thinned <- thin(loc.data = spoccs, lat.col = "y", long.col = "x", spec.col = "Species", 
    thin.par = 5, reps = 1000, verbose = T, out.dir = "./Thinned", locs.thinned.list.return = TRUE, 
    write.files = TRUE, out.base = "Petromarula pinnata thinned")
## ********************************************** 
##  Beginning Spatial Thinning.
##  Script Started at: Sun May 17 22:17:19 2020
## lat.long.thin.count
##   39 
## 1000 
## [1] "Maximum number of records after thinning: 39"
## [1] "Number of data.frames with max records: 1000"
## [1] "Writing new *.csv files"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin1.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin2.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin3.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin4.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin5.csv"

6 Μοντελοποίηση

Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο εδώ

6.1 Εισαγωγή

Είμαστε πλεόν στην ευχάριστη θέση να έχουμε στην διάθεση μας όλα τα δεδομένα τα οποία είναι απαραίτητα για να προχωρήσουμε στο κυρίως σκέλος της ανάλυσης μας. Κανονικά θα έπρεπε να χρησιμοποιήσουμε πολλά μοντέλα, τα οποία ανήκουν στην ίδια ομάδα (Παρουσίας/Απουσίας ή Ψεύδο-απουσίας - Presence/Absence: P/A), αλλά και σε άλλη ομάδα (Μόνο Παρουσίας - Presence Only: P/O). Κάθε ένα εκ των μοντέλων αυτών έχει τα πλεονεκτήματα και τα μειονεκτήματα του14. Στο σημείο αυτό, πρέπει να τονίσουμε ότι είναι εντελώς λάθος να χρησιμοποιούμε δύο μοντέλα15 τα οποία ανήκουν σε διαφορετικές ομάδες (εν αντιθέσει με ότι προτείνεται στο παράδειγμα του προαναφερθέντος βιβλίου), αλλά να χρησιμοποιεί πολλά μοντέλα τα οποία ανήκουν στην ίδια ομάδα (π.χ., P/A), προκειμένου να μπορεί να βγάλει ασφαλή συμπεράσματα σχετικά με το ερώτημα το οποίο καλείται να απάντησει και να είναι σίγουρος ότι δεν πρόκειται περί στατιστικών κατασκευασμάτων (Araújo & New, 2007).

6.2 Φόρτωση των δεδομένων θέσης

Τώρα θα φορτώσουμε ένα από τα set δεδομένων16 τα οποία έχουν προέλθει από το thinning των θέσεων εμφάνισης17.

thinned_1 <- readxl::read_excel("./Thinned/Petromarula pinnata thinned_thin1.xlsx")

6.3 Μορφοποιήση των δεδομένων

Το πρώτο πράγμα το οποίο θα χρειαστεί να κάνουμε, είναι να μορφοποιήσουμε με τον κατάλληλο τρόπο τα δεδομένα μας, ώστε να μπορέσουμε να τρέξουμε τις εντολές από τα πακέτα biomod και ecospat. Το δεύτερο πακέτο θα το χρησιμοποιήσουμε στην περίπτωση όπου οι θέσεις εμφάνισης του είδους που μας ενδιαφέρει είναι \(<50\)18.

Ένα άλλο σημείο το οποιο πρέπει να προσέξουμε, έχει να κάνει με τον αριθμό των απουσιών (absences - στην πλειονότητα των περιπτώσεων βέβαια πρόκειται για ψευδοαπουσίες -pseudoabsences, καθώς στην πραγματικότητα δεν έχουμε δεδομένα απουσίας του εκάστοτε είδους). Ο αριθμός αυτός πρέπει να είναι ίσος με τον αριθμό των παρουσιών19.

Τέλος, ένα τρίτο, κρίσιμο σημείο είναι ο αριθμός των set των pseudoabsences που θα περιλαμβάνει το μοντέλο μας, προκειμένου να αποφύγουμε τον σκόπελο της μεροληπτικής δειγματοληψίας (sampling bias)20.

sp_formated_data <- BIOMOD_FormatingData(resp.var = rep(1, nrow( thinned_1 ) ),
                                         
                                    ## Our rasterStack with the uncorrelated
                                    ## variables
                                    expl.var = current_climate_Crete,
                                    
                                    ## The names of our coords
                                    resp.xy = thinned_1[,c('x', 'y')],
                                    
                                    ## Species' name
                                    resp.name = "Petromarula pinnata",
                                    
                                    ## Number of pseudoabsences' sets
                                    PA.nb.rep = 2,
                                    
                                    ## Number equaling that of the presences
                                    PA.nb.absences = nrow(thinned_1), 
                                    
                                    ## You can change this to 'disk' if you'd
                                    ## like
                                    PA.strategy = 'random',
                                    
                                    ## Min distance from occurrences
                                    PA.dist.min  =  NULL,
                                    
                                    ## Max distance from occurrences
                                    PA.dist.max  =  NULL) 
## 
## -=-=-=-=-=-=-=-=-=-= Petromarula pinnata Data Formating -=-=-=-=-=-=-=-=-=-=
## 
##  Response variable name was converted into Petromarula.pinnata
##       ! No data has been set aside for modeling evaluation
##    > Pseudo Absences Selection checkings...
##    > random pseudo absences selection
##    > Pseudo absences are selected in explanatory variables
##          ! Some NAs have been automaticly removed from your data
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
sp_formated_data
## 
## -=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.formated.data.PA' -=-=-=-=-=-=-=-=-=-=-=-=
## 
## sp.name =  Petromarula.pinnata
## 
##   31 presences,  0 true absences and  76 undifined points in dataset
## 
## 
##   4 explanatory variables
## 
##      bio_13          bio_15          bio_3           bio_5      
##  Min.   :106.0   Min.   :78.00   Min.   :29.00   Min.   :209.0  
##  1st Qu.:138.0   1st Qu.:83.00   1st Qu.:30.00   1st Qu.:269.0  
##  Median :153.0   Median :84.00   Median :31.00   Median :277.0  
##  Mean   :159.6   Mean   :83.82   Mean   :31.08   Mean   :274.3  
##  3rd Qu.:187.0   3rd Qu.:85.00   3rd Qu.:32.00   3rd Qu.:284.0  
##  Max.   :213.0   Max.   :88.00   Max.   :33.00   Max.   :301.0  
## 
## 
##  2 Pseudo Absences dataset available ( PA1 PA2 ) with  39 
## absences in each (true abs + pseudo abs)
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

Ας αναπαραστήσουμε γραφικά ότι φτιάξαμε στο προηγούμενο βήμα:

plot(sp_formated_data)

Η εικόνα αυτή μας δείχνει που βρίσκονται οι ψευδοαπουσίες (pseudoabsences) στο set των ψευδοαπουσιών σε σχέση με τις θέσεις εμφάνισης. Όπως αναμενόταν, οι ψευδοαπουσίες είναι τυχαία κατανεμημένες.

6.4 Παραμετροποιήση των μοντέλων

Πλέον μπορούμε να παραμετροποιήσουμε τα μοντέλα που θέλουμε να τρέξουμε:

modelling_options <- BIOMOD_ModelingOptions()

Τώρα, προσδιορίζουμε ποιά μοντέλα θα τρέξουμε, καθώς και πολλές άλλες λεπτομέρειες:

## Always change the name of this object, so as to match the name of your species
Petromarula_pinnata_models <- BIOMOD_Modeling( data = sp_formated_data,
                                       
                                       
                                       ## Here, we indicate which models we are
                                       ## going to run (remember that they need
                                       ## to be from the same group P/A ή P/O).
                                       ## Here, we are going to use some models
                                       ## belonging to the P/A group
                                       ## ('GLM','GBM','GAM','CTA','ANN','SRE',
                                       ## 'FDA','MARS','RF')
                                       models  =  c('RF', 'CTA'),
                                       
                                       models.options = modelling_options,
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate its'
                                       ## results [min: 10 - usually: 100]
                                       NbRunEval = 1,
                                       
                                       ## Training-Testing data split
                                       DataSplit = 80, 
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate the
                                       ## importance of each IV [min: 10]
                                       VarImport = 1,
                                       
                                       ## Evaluation criteria
                                       models.eval.meth  =  c('TSS','ROC','KAPPA'), 
                                       
                                       ## Whether or not we will save the results
                                       SaveObj  =  TRUE,
                                       
                                       ## Should the models be re-scaled in 0-1000 range?
                                       rescal.all.models  =  TRUE,
                                       
                                       ## Should our models be weighted based on all the data?
                                       do.full.models  =  TRUE, 
                                       
                                       ## The name of our analysis
                                       modeling.id  =  "allmodels")

Ας φτιάξουμε τους φακέλους όπου θα αποθηκεύσουμε τα μοντέλα μας:

## Create some folders to store our results
dir.create("./BIOMOD2 Modelling Outpout/", recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/", recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata", recursive = TRUE, 
    showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current", 
    recursive = TRUE, showWarnings = FALSE)


## Save the results in rds format
saveRDS(Petromarula_pinnata_models, "./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current/ Petromarula_pinnata models.rds")

6.5 Αξιολόγηση των αρχικών μοντέλων

Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων:

Petromarula_pinnata_scores <- get_evaluations(Petromarula_pinnata_models)

6.6 Σημαντικότητα μεταβλητών

Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA), κάθε σετ δεδομένων (PA1-2) και κάθε επανάληψη (Run, Full):

(Petromarula_pinnata_var_import <- get_variables_importance(Petromarula_pinnata_models))

Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA) συνολικά (δηλαδή τον μέσο όρο):

options(digits = 2)
apply(Petromarula_pinnata_var_import, c(1, 2), mean)
##           RF   CTA
## bio_13 0.528 0.995
## bio_15 0.091 0.069
## bio_3  0.093 0.088
## bio_5  0.344 0.373

6.7 Ensemble modelling

Τώρα θα τρέξουμε την ανάλυση που θα ενώνει τα αποτελέσματα των δύο μοντέλων (ονομάζεται ensemble modelling αυτή η τεχνική ανάλυσης):

Petromarula_pinnata_ensemble_models <- BIOMOD_EnsembleModeling( modeling.output = Petromarula_pinnata_models,
                                                   em.by = 'all',
                                                   chosen.models  =  'all',
                                                   
                                                   # You can change that
                                                   eval.metric = 'TSS', 
                                                   
                                                   # Keep only the models with
                                                   # TSS >= 0.8
                                                   eval.metric.quality.threshold = 0.8, 
                                                   models.eval.meth = c('KAPPA','TSS','ROC'),
                                                   prob.mean  =  T, 
                                                   committee.averaging = TRUE,
                                                   prob.mean.weight  =  T, 
                                                   prob.mean.weight.decay  =  'proportional',
                                                   VarImport = 0 )
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_ensemble_models, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble/Petromarula_pinnata ensemble.rds")

6.8 Αξιολόγηση των συνολικών μοντέλων (ensemble modelling)

Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων από το ensemble modelling. Όσο πιο κοντά στη μονάδα είναι οι τιμές για τους δείκτες TSS, ROC και KAPPA, τόσο καλύτερα21:

options(digits = 4)
(Petromarula_pinnata_ensemble_models_scores <- get_evaluations(Petromarula_pinnata_ensemble_models))
## $Petromarula.pinnata_EMmeanByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.864  636.0       90.32       96.05
## TSS          0.864  636.0       90.32       96.05
## ROC          0.980  637.5       90.32       96.05
## 
## $Petromarula.pinnata_EMcaByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.861    747        87.1       97.37
## TSS          0.845    747        87.1       97.37
## ROC          0.958    750        87.1       97.37
## 
## $Petromarula.pinnata_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.864    641       90.32       96.05
## TSS          0.864    641       90.32       96.05
## ROC          0.982    640       90.32       96.05
## ===========================================================##
## Vizualisation
## ===========================================================##

g2 <- models_scores_graph(Petromarula_pinnata_ensemble_models, by = "models", 
    metrics = c("TSS", "KAPPA"))

6.9 Προβολή στο μέλλον

Ας τρέξουμε την ανάλυση για το πώς θα αλλάξει η κατάσταση του είδους που μας ενδιαφέρει στο μέλλον. Αρχικά θα τρέξουμε την ανάλυση για κάθε μοντέλο ξεχωριστά (BIOMOD_Projection) και συνολικά (BIOMOD_EnsembleForecasting) για τις παρούσες συνθήκες και στη συνέχεια για τις μελλοντικές συνθήκες (με τις ίδιες εντολές, αλλά θα αλλάξουμε το raster αρχείο στην εντολή new.env):

Petromarula_pinnata_models_proj_current <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models, 
    new.env = current_climate_Crete, proj.name = "current", binary.meth = "TSS", 
    selected.models = "all", build.clamping.mask = TRUE, output.format = ".img", 
    do.stack = T)
Petromarula_pinnata_ensemble_models_proj_current <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models, 
    projection.output = Petromarula_pinnata_models_proj_current, binary.meth = "TSS", 
    output.format = ".img", do.stack = T)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual/Petromarula_pinnata Current Individual Projections.rds")

saveRDS(Petromarula_pinnata_ensemble_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble/Petromarula_pinnata Current Ensemble Projections.rds")
Petromarula_pinnata_proj_2070_cc85 <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models, 
    new.env = future_climate_Crete, proj.name = "2070_CC85", binary.meth = "TSS", 
    output.format = ".img", do.stack = T)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070", recursive = TRUE, 
    showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual/Petromarula_pinnata 2070 Individual Projections.rds")
Petromarula_pinnata_ensemble_models_proj_2070_cc85 <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models, 
    projection.output = Petromarula_pinnata_proj_2070_cc85, binary.meth = "TSS", 
    output.format = ".img", do.stack = T)
saveRDS(Petromarula_pinnata_ensemble_models_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds")

6.10 Οπτικοποιήση των συνολικών μοντέλων

Τέλος, ας οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας. Πρώτα θα θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει την αξιολόγηση της κατάστασης της Petromarula pinnata μέσω της τεχνικής ensemble modelling. Στη συνέχεια, θα κρατήσουμε μόνο τα raster layers για το mean και το weighted mean. Τέλος, θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει μόνο τις συντεταγμένες των θέσεων εμφάνισης του είδους που μας ενδιαφέρει:

stk_Petromarula_pinnata_ef_2070_cc85 <- get_predictions(Petromarula_pinnata_ensemble_models_proj_2070_cc85)

stk_Petromarula_pinnata_ef_2070_cc85 <- subset(stk_Petromarula_pinnata_ef_2070_cc85, 
    grep("EMmean|EMwmean", names(stk_Petromarula_pinnata_ef_2070_cc85)))

names(stk_Petromarula_pinnata_ef_2070_cc85) <- sapply(strsplit(names(stk_Petromarula_pinnata_ef_2070_cc85), 
    "_"), getElement, 2)

Petromarula_pinnata_coords <- thinned_1 %>% dplyr::select(x, y)

Στη συνέχεια θα οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας έως αυτό το σημείο:

library(viridis)
levelplot(stk_Petromarula_pinnata_ef_2070_cc85, main = "Petromarula pinnata ensemble projections\nin 2070 with cc85", 
    col.regions = viridis(256)) + layer(sp.points(SpatialPoints(Petromarula_pinnata_coords), 
    pch = 20, col = 1))

6.11 Μεταβολή περιοχής εξάπλωσης

Έχουμε φτάσει πλέον στο τελικό σημείο της ανάλυσης μας, όπου θα φορτώσουμε τα αποτελέσματα της ανάλυσης για το mean και το weighted mean της Petromarula pinnata και στη συνέχεια θα υπολογίσουμε και θα οπτικοποιήσουμε τη μεταβολή της περιοχής εξάπλωσης της στο μέλλον:

##==================================
## Load the the data
##==================================

## Current
Petromarula_pinnata_bin_proj_current <- stack('Petromarula.pinnata/proj_current/proj_current_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_current) <- c('ca', 'm', 'wm')

## 2070
Petromarula_pinnata_bin_proj_2070_CC85 <- stack('Petromarula.pinnata/proj_2070_CC85/proj_2070_CC85_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_2070_CC85) <- c('ca', 'm', 'wm')

##=================================
## Species Range Change
##=================================

SRC_current_2070_CC85 <- BIOMOD_RangeSize( Petromarula_pinnata_bin_proj_current,
                                           Petromarula_pinnata_bin_proj_2070_CC85 )
SRC_current_2070_CC85$Compt.By.Models
##    Loss Stable0 Stable1 Gain PercLoss PercGain SpeciesRangeChange
## ca  106     250      25   85   80.916   64.885            -16.031
## m    95     300      14   57   87.156   52.294            -34.862
## wm  110     250      23   83   82.707   62.406            -20.301
##    CurrentRangeSize FutureRangeSize.NoDisp FutureRangeSize.FullDisp
## ca              131                     25                      110
## m               109                     14                       71
## wm              133                     23                      106
Petromarula_pinnata_src_map <- SRC_current_2070_CC85$Diff.By.Pixel
names(Petromarula_pinnata_src_map) <- c("ca cur-2070", "m cur-2070", "wm cur-2070")

##=================================
## Vizualisation
##=================================

my.at <- seq(-2.5,1.5,1)

myColorkey <- list(
                   ## color change
                   at = my.at, 
                   labels = list(
                     
                     ## labels
                     labels = c("lost", "pres", "abs","gain"),
                     
                     ## legend
                     at = my.at[-1]-0.5 
                   ))

rasterVis::levelplot( SRC_current_2070_CC85$Diff.By.Pixel,
                      main = "Petromarula pinnata range change",
                      colorkey = myColorkey,
                      layout = c(1,3),
                      )

7 Εργασία για το σπίτι

Έχετε διορία δεκατεσσάρων (14) ημερών να στείλετε την εργασία σας σε μορφή PDF22 σε αυτό το e-mail.
1. Κατεβάστε γεωγραφικά δεδομένα για την Ελλάδα.
2. Δημιουργήστε ξεχωριστά αρχεία για την Κρήτη, την Αττική, τον Άθω, την Μακεδονία και τη Θεσσαλία.
3. Κατεβάστε τα απαραίτητα κλιματικά δεδομένα.
4. Κατεβάστε υψομετρικά δεδομένα για την Ελλάδα.
5. Υπολογίστε τις βασικές στατιστικές παραμέτρους για όλα τα αρχεία τα οποία δημιουργήσατε.
6. Υπολογίστε τον συντελεστή συσχέτισης για όλα τα περιβαλλοντικά δεδομένα για όλα τα αρχεία που έχετε δημιουργήσει. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.7\), θα τις χρησιμοποιήσετε όλες;
7. Κατεβάστε δεδομένα θέσης για τα είδη Asperula pubescens, Centaurea idaea και Scutellaria sieberi.
8. Παραθέστε πληροφορίες σχετικά με τα είδη αυτά (Καθεστώς κινδύνου/προστασίας, καθεστώς ενδημισμού, ταξινομικές πληροφορίες). Προχωρήστε στις απαραίτητες ενέργειες, έτσι ώστε να διαπιστώσετε εάν τα δεδομένα αυτά εντοπίζονται στη σωστή γεωγραφικά περιοχή.
9. Πόσες θέσεις εμφάνισης αριθμούσαν τα είδη αυτά πριν και μετά την διαδικασία αραίωσης;
10. Φορτώστε τα τελικά δεδομένα θέσης για κάθε είδος για το οποίο έχετε κατεβάσει δεδομένα θέσης.
11. Αλλάξτε τα rasterStack αρχεία με αυτά για την Κρήτη σε όλα τα απαραίτητα σημεία του κώδικα. Αλλάξτε επίσης το όνομα του είδους σας. Αλλάξτε τον αριθμό των σετ σε 10 και τον αριθμό των ψευδοαπουσιών23. Αλλάξτε τον αριθμο των φορών που θα τρέξει η ανάλυση σε 10.
12. Ποιά είναι η σημαντικότερη μεταβλητή σε κάθε αλγόριθμο; Συμπίπτουν μεταξύ των αλγορίθμων;
13. Ποιό είναι το καλύτερο μοντέλο σύμφωνα με τους δείκτες αξιολόγησης TSS και KAPPA;.
14. Προβάλετε τα μοντέλα σας στο μέλλον. Παραθέστε τους σχετικούς χάρτες.
15. Πώς μεταβάλλεται η περιοχή εξάπλωσης των ειδών αυτών στο μέλλον;24 Τι συμβαίνει στην περίπτωση της μηδενικής και της πλήρους διασποράς, αντίστοιχα;
16. Παραθέστε το σχετικό γράφημα.

Βιβλιογραφία

Araújo, M.B. & New, M. (2007) Ensemble forecasting of species distributions. Trends Ecol. Evol., 22, 42–47.

Guisan, A., Thuiller, W., & Zimmermann, N.E. (2017) Habitat Suitability and Distribution Models with Applications in R. Cambridge University Press,

McSweeney, C.F., Jones, R.G., Lee, R.W., & Rowell, D.P. (2015) Selecting CMIP5 GCMs for downscaling over multiple regions. Clim. Dyn., 44, 3237–3260.

  1. Υπάρχουν δύο τρόποι για να το κάνετε αυτό γρήγορα:
    1. library(easypackages) packages('package_name_1', 'package_name_2')
    2. install.packages(c('package_name_1', 'package_name_2'), dependecies = TRUE)↩

  2. Ή με αυτήν την εντολή εάν θέλουμε μόνο την Ελλάδα: ISO_countries <- getData("ISO3") %>% as.data.frame %>% filter(NAME=="Greece")↩

  3. Ένα τυπικό παράδειγμα είναι ο Άτλαντας της Χλωρίδας του Αιγαίου ή ο Άτλαντας της Χλωρίδας της Ευρώπης, όμως στην περίπτωση αυτή υπεισέρχονται άλλα προβλήματα τα οποία σχετίζονται με την ανάλυση των εν λόγω χαρτών (π.χ., ο Άτλαντας της Χλωρίδας της Ευρώπης είναι σε ανάλυση \(50 \times 50 km^2\))↩

  4. Δεν είναι τα μοναδικά όμως: υπάρχουν πολλά άλλα, όπως το dismo και το spocc↩

  5. Κάθε ένα από αυτά τα πακέτα, έχει τουλάχιστον μια σύντομη περιγραφή (αγγλιστί: vignette), την οποία είναι ΑΠΑΡΑΙΤΗΤΟ να διαβάσετε - ειδικά αυτές του rgbif↩

  6. Με την εντολή ?rgbif::name_suggest() μπορείτε να δείτε τα arguments που δέχεται η εντολή αυτή και τι μπορείτε να αλλάξετε - σας συνιστώ εντόνως να το κάνετε, αλλά και για κάθε εντολή που χρησιμοποιείτε↩

  7. Καλό είναι να τις διαβάσετε, γιατί θα σας βοηθήσει και στην εργασία σας↩

  8. Και εγώ στη θέση σας, το ίδιο θα έκανα↩

  9. Γιατί;↩

  10. Αυτό είναι το λεγόμενο hard boundary - άλλοι ερευνητές εφαρμόζουν το όριο \(> 0.8\) ή ακόμα και \(> 0.9\), αλλά είναι επιστημονικά ασφαλέστερο να εφαρμόζει κανείς το αυστηρότερο όριο↩

  11. Εδώ εφαρμόζουμε το αυστηρότερο όριο - μπορείτε να το αλλάξετε, εάν επιθυμείτε↩

  12. Σε απόσταση \(<1000m\) όταν εργαζόμαστε με χάρτες υποβάθρου κλίμακας 1 x 1 km2 (ή ανάλυσης 30 arc-seconds)↩

  13. Θα σας βοηθήσει το vignette του εν λόγου πακέτου↩

  14. Προτρέξτε στους Guisan et al. (2017) για περισσότερες λεπτομέρειες - εάν χρειάζεστε περισσότερη βιβλιογραφία, απευθυνθείτε σε εμένα (εάν φυσικά το επιθυμείτε)↩

  15. Όπως θα διαπιστώσετε παρακάτω, θα εφαρμόσουμε μια τεχνική η οποία ονομάζεται ensemble modelling. Σύμφωνα με αυτή, είναι προτιμότερο και επιστημονικά ορθότερο να μην χρησιμοποιεί κανείς ένα μόνο μοντέλο (π.χ., GLM ή τo MaxEnt)↩

  16. Προσέξτε ότι πρόκειται για αρχείο xlsx και όχι csv - αυτό σημαίνει ότι άνοιξα στο excel το αρχείο csv, χρησιμοποίησα την εντολή data to columns και το έσωσα ως αρχείο xlsx↩

  17. Δεν έχει σημασία ποιό θα διαλέξουμε - έχουν όλα την ίδια πιθανότητα↩

  18. Θα μπορούσε κάποιος να αντιτείνει ότι μπορούμε να χρησιμοποιήσουμε το πακέτο ecospat όταν \[\frac{\sum{Ανεξάρτητες,\, μη\, συγγραμικές\, μεταβλητές}}{\sum{Θέσεις\, εμφάνισης}} \geq 0.1\]↩

  19. Στο βιβλίο των Guisan et al. (2017) αναφέρεται λανθασμένα ότι ο αριθμός των pseudoabsences πρέπει να είναι τουλάχιστον 10000 - ο αριθμός αυτός έχει να κάνει με το MaxEnt, το οποίο ανήκει στην ομάδα των P/O μοντέλων↩

  20. Καλό είναι να είναι τουλάχιστον 10 τα set αυτά↩

  21. Θυμηθείτε να με ρωτήσετε τι σημαίνουν οι υπόλοιπες μεταβλητές↩

  22. Διαφορετικά δεν θα γίνει δεκτή και δεν θα βαθμολογηθείτε↩

  23. Πάντα σύμφωνα με το είδος σας↩

  24. Παραθέστε ΟΛΑ τα στοιχεία από τον σχετικό πίνακα↩

Κώστας Κουγιουμουτζής

Spring semester 2019-2020

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr